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Abstract 


A fourth-order and rapid numerical algorithm, utilizing a procedure as 
Runge-Kutta methods, is derived for solving nonlinear equations. The 
method proposed in this article has the advantage that it, requiring no 
calculation of higher derivatives, is faster than the other methods with the 
same order of convergence. The numerical results obtained using the devel- 
oped approach are compared to those obtained using some existing iterative 
methods, and they demonstrate the efficiency of the present approach. 
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1 Introduction 


The development of numerical techniques for solving nonlinear (algebraic 
or transcendental) equations (NEs) is a subject of considerable interest. A 
vast amount of research work has been invested in the study of the NEs. 
Solution of these equations can be obtained using classical numerical methods 
as Newton—Raphson method (NRM) or the Householder iteration method 
(HIM); see [6]. There are many articles that deal with NEs, for example, 
[1, 2, 3, 4, 5]. A more extensive list of references as well as a survey on 
the progress made on this class of problems may be found in [8]. However, 
in order to establish a more accurate algorithm, one has to calculate higher 
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order derivatives, for example, like the HIM, which can be considered a serious 
drawback. 

Here, we establish a fourth-order fast algorithm for finding an accurate 
solution of the NEs with the advantage that there is no calculation of higher 
derivatives. Moreover, the merit in this method can be observed due to high 
exactness, lesser number of iterations, and lesser value of errors as compared 
to the other existing methods. The examples analyzed here reveal that the 
developed algorithm is effective to solve the NEs. 


2 Basic idea of the new algorithm 


Consider the nonlinear algebraic equation 


f(z) =0. (1) 


We assume that a is a simple zero of (1), that f is a C? function on an 
interval containing a, and that A is an initial guess sufficiently close to a. 
Also we suppose that |f’(a)| > 0. Using Taylor’s series around A for (1), we 
have 

1 1 


POA) + @ AFA) + 5 (@— APIA) + G@—APIMOA) +2 = 0. (2) 


We can rewrite (2) in the following form: 


x=c+N(a), (3) 
where FO) 
=— ay (4) 
and F aN 
N(e) =-5(e PE (5) 
The NRM is then given by 
ae (@n) 
Int1 = In fl(@n)’ (6) 


where 2%o is the initial guess. The iteration (6) will converge to a if the 
starting point x9 is close enough to a. This process has the second-order 
convergence. 

The HIM is given by (sce [6]) 


f (ix) i (tn) f" (fn) 
Bett =n Fi) OF) my 
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The iteration (7) with the third-order convergence will converge to a if 
the starting point xp is close enough to a. 


As observed in the above-mentioned methods, for obtaining a more ac- 
curate solution of a, we have to take n, order iteration, as large as possible, 
or go forward with the calculation of higher derivatives. In order to remove 
these demerits, here we look for a scheme without computing higher order 
iterations/derivatives that contain high accuracy. To illustrate the idea be- 
hind the developing algorithm, we construct an iterative procedure in the 
following form: 


Tn+1 =I + Ak, + Bko, (8) 
where fen) 
In 
ky = fi(tp)’ (9) 
_ flan +Ch) 
kp = Pan) (10) 


In (8)—(10), the coefficients A, B, and C are constants to be determined 
later. We note that the NRM (6) can be obtained from (8) for the values 
A=-land B=C=0. 

Now, we want to determine the constants of A, B, and C so that the 
formula (8) is similar to the formula (7) of accuracy. From (8), therefore, 
using Taylor’s expansion, we have 


Inti = On +(A+B4 po) fee) - Bort ee (11) 


Comparing (7) and (11) and equating the coefficients with identical terms, 

we get 
A+B+BC=-1, 
{nee =-1. 2) 
The solutions of (12) allow us to suggest different iterative methods (of 
at least third-order) for solving the NEs of (1). Two solutions for the system 
(12) are A=0, B= atv6 C= 1 V5 and A= B=C=~-1. Thus, sub- 
stituting the later solution into (8)—(10), we will have the following iterative 

algorithm: 


— f(@n) 
ae C8 (13) 
_ f (an = ky) 
ko ; 
Now, if, according to the definition f’(x) = limp+o fet hy He) we esti- 


mate the value of first derivative by placing h ~ — a as 
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! Ln tn) — In —- (an) 
wey w Lol =a res) as 


then we will have 
Lnt1 = Ln — [ki + ke], 
f?(@n) 
F'(en) (Flen) — f (an - 3) ae 


Now the second algorithm could be gained by substituting the former 
solution into (8)—(10) as 


3475 
Int1 = Ty — 3a Kas 
ky = 
1 fl(@n)’ (16) 
k f(&n =F 1—V5 ,) 
_ f' (zn) , 


Remark 1. As we know well, these kinds of algorithms need a good initial 
guess to work. This shortcoming can be readily removed by finding subinter- 
vals on x that contain sign changes of f(x) (observe the appendix section). 


Remark 2. From (15) and (16), we can comprehend that in each iteration 
of the algorithm (15), we need to perform four function evaluations, and for 
the algorithm (16), there are three function evaluations. However, the fast 
convergence of these methods often compensates for these defects. 


3 Analysis of convergence 


Before proving the convergence of the algorithms (15) and (16), here the 
following definition is given. 


Definition 1. Let e, = x, —a be the truncation error in the nth iterate. If 
there exist a number 8 > 1 and a constant c 4 0 such that 


lim lens] =c (17) 


noo |e,|F i 


then ( is called the order of convergence of the method. 


Now, we consider the convergence of the algorithms (15) and (16). 
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Theorem 1. Consider the nonlinear algebraic equation f(x) = 0. Suppose 
that f is sufficiently differentiable. Then the convergence of the algorithm 
(16) is at least of order 3. 


Proof. Let a be a simple zero of f(a). Since f(a) is sufficiently differentiable, 
by expanding f(z), f’(x), f’ (a), and f’’(a) around a, we get 


f(an) = f'(@) [en + doe? + dze3 + dae + doen +---], 

f'(&n) = f'(a) [1 + 2dzen + Bdge? + 4dase3 + 5d5e%, + 6deeg +--+], 

f'" (an) = f' (a) [2dz + 6dgen + 12dye? + 20dse3 + 30dge4 + 42dze> +--+], 

f'" (tn) = f' (a) [6d3 + 24daen + 60d5e? + 120dge3 + 210d7e4 + 336dge> +--+], 


(18) 
where d,, = 4 ao} n = 2,3,..., and e) =X», — a. We can rewrite (16) as 
follows: Han) 

Ln 
n = Cn B n ’ 
Cn+1 = en + aG +offn) (19) 


or using Taylor’s series as 


a+ eysten)+ F (Le ) Pet = (fe ren) | 


E€n+1 = En +B 


2 J’ ) 6 F' ay,) 

(20) 
where B = sive and C = By a simple operation, from (18) and 
(20), we acquire 

7 6 fl2 5 — 3)f! mr 


6 f(a) 


which shows that the algorithm (16) is at least a third-order convergent 
method and this ends the proof. 


Theorem 2. Under the assumptions of Theorem 1, the convergence of the 
algorithm (15) is at least of order 4. 


Proof. Consider the above expansions of f(x), f’(x), f(x), and f’’(x). Pro- 
ceeding as before, we obtain 


ax leagil 23 3 f(a) + f(a) f (a) 
naeo ent? 24 f(a) 


(22) 


which shows that the algorithm (15) is at least a fourth-order convergent 
method and thus the theorem is proved. 
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4 Numerical implementations 


The test problems given in this section demonstrate the effectiveness of the 
developed algorithms to solve the NEs. Meanwhile, the results obtained using 
the algorithms proposed in this work, (15) and (16) are compared to those 
obtained using the HIM and a fourth-order method presented in [7]. All 
calculations have been carried out by using the symbolic calculus software 
Maple 17 to 2500 significant decimal digits. 


Example 1. Solve the following NE: 


x — cos(x) = 0, Lo = 2. (23) 


Example 2. Solve the following NE: 


z-2-e *=0, Lo = 2. (24) 


Example 3. Solve the following NE: 


sin? (x) — 27 +1 =0, Xo = —2. (25) 


Example 4. Solve the following NE: 


=a)" SU, to =. (26) 


The numerical results of Examples 1—4 can be observed in Table 1. In 
Table 1, we list the value of absolute error f(x») (labeled as |f(xn)|) when 
n = 1,2,3,4,5 for the above examples. 


From the numerical results of these cases analyzed here, it is easy to 
conclude that the proposed algorithms are effective and accurate for solving 
NEs. Moreover, unlike HIM, the developed approaches require just first-order 
derivative per step. 


5 Conclusions 


The presented algorithms in this article give rapid convergent approxima- 
tions. They have the advantage of giving an accurate solution with fewer 
iterations than the other existing methods. 
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Table 1: The numerical results obtained from solving the above investigated examples 
using the algorithms (15) and (16), HIM, and the method of [7] for few iterations. 


Algorithm lf(zi1)|_ | F (x2) lf(zs)|_|f(ea)l__ | f(s) | 
Example 1 

HIM 7.6_3 7.7_¢6 7.812 8.024 8.5_48 

Algorithm (16) 1.14 5.0_5 5.6_15 7.7_45 2.0_134 


Method of [7] 9.9_4 1.244 2.4_58 4.1_933 3.3932 
Algorithm (15) 1.24 9.1_19 1.6_76 1.4396 9.4—1227 


Example 2 
HIM G04 Ad_e Bilge Siow ABce 
Algorithm (16) 1.3_6 1.9_91 5.8_66 1.7_199 4.2600 
Method of [7] 7.9_8 9.3_33 1.9132 2.9531 1.6_2126 
Algorithm (15) 4.1_8 3.234 1.1_138 1.5_556 5.5_2298 


Example 3 
HIM 38-4 R26 80.4  O0se “27x 
Algorithm (16) 1.6_1 5.6_4 3.311 7.0_33 6.7_98 
Method of [7] 6.9_2 2.2_6 2.894 7.3_96 3.0382 
Algorithm (15) 3.6_2 4.8_8 1.7_31 2.6_125 1.3_500 
Example 4 
HIM 09.4 Ot 5 304° Gl»  -2.6-%6 
Algorithm (16) 5.1_2 6.0_5 1.3_13 1.1_39 8.7_118 
Method of [7] 1.6_2 4.42 2.4_ 39 2.0_119 1.0_475 


Algorithm (15) 2.9_4 2.315 9.9_60 3.3237 4.1_947 


Appendix 


Maple code of the algorithm (15) for finding a simple zero of f(x) := fun with 
the initial guess xo := x0 = (a,b) up to dgt digits, where ftol determines 
the tolerance. 


fRoots := proc(fun,x0,ftol,dgt) 
local m, n, p, r, i, j, ki, k2, f, Df, df, idx, num, fnrm, it; 
if dgt<=abs(log10(ftol)) then 
error " Number of digits must be greater than log of ftol" 
end if; 
idx := [°$°(x0[1]+ii/2,ii = 0 .. 2*(x0(2]-x0[1]))]; 
f := unapply(fun,x); 
Df := diff (f(x),x); 
df := unapply(Df,x); 
num := 2*(x0[2]-x0[1]); 
m := 0; 
for i to num-1 do 
if signum(f(idx[i])) <> signum(f(idx[it1])) then 
m := m+1; 
pl[m) := [idx[i], idx[i+1]]; 
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end if 
end do; 
if not type(p,table) then error " f(x) contains no sign changes " end if; 
for j from 1 to m do 
r({0O] := 1/2*(p[j] [1] +plj1[2]); 
fnrm : =1: it := 0: 
for n from 0 while (fnrm > ftol) do 
k1 := evalf(f(r[n])72/(df (r[n] )* (f(r [n] )-f (x [n] -f (x [n]) /df (r[n])))) ,dgt) ; 
k2 := evalf(ki*f(r[n]-k1)/£(r[n]) ,dgt) ; 
r[n+1] := evalf(r[n]-k1-k2,dgt) ; 
fnrm := abs(evalf((f(r[nt1])) ,dgt)); 
it := it+1; 
end do; 
printf(" x[0] = %g \n Itration = 4d \n |f£(x)| = he 
x = Za \n\n",r[0],it,fnrm,r[it]); 
end do; 
end proc: 


For example, one can use the above code to solve the NE 1n(x+1)+x-1 
(with root « = 0.557145598997611416858672000001) by the following: 


fRoots(1n(xt+1)+x-1, [0,10] ,107(-10) ,30); 
x[0] = 0.75 
Itration = 2 
l£(x) | = 6.839800e-26 

x = .557145598997611416858671958351 
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